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We study the impact of the Peterlin approximation on the statistics of the end-to-end separation 
of polymers in a turbulent flow. The FENE and FENE-P models are numerically integrated along 
a large number of Lagrangian trajectories resulting from a direct numerical simulation of three- 
dimensional homogeneous isotropic turbulence. Although the FENE-P model yields results in qual¬ 
itative agreement with those of the FENE model, quantitative differences emerge. The steady-state 
probability of large extensions is overestimated by the FENE-P model. The alignment of polymers 
with the eigenvectors of the rate-of-strain tensor and with the direction of vorticity is weaker when 
the Peterlin approximation is used. At large Weissenberg numbers, both the correlation times of 
the extension and of the orientation of polymers are underestimated by the FENE-P model. 


I. INTRODUCTION 

The addition of elastic polymers to a Newtonian solvent introduces a history dependence in the response of 
the fluid to a deformation and hence modifies the rheological properties of the solvent A In turbulent flows, the 
non-Newtonian nature of polymer solutions manifests itself through a considerable reduction of the turbulent 
drag compared to that of the solvent alone^— What renders this phenomenon even more remarkable is that 
an appreciable drag reduction can already be observed at very small polymer concentrations (of the order of 
a few parts per million). Turbulent drag reduction was discovered by Tomsi more than sixty years ago and 
is nowadays routinely used to reduce energy losses in crude-oil pipelines^ A full understanding of turbulent 
drag reduction nonetheless remains a difficult challenge, because in the turbulent flow of a polymer solution 
the extensional dynamics of a large number of polymers is coupled with strongly nonlinear transfers of 
kinetic energy. 

The study of turbulent flows of polymer solutions is essentially based on two approaches: the molecular 
approach and the continuum one. In the molecular (or Brownian Dynamics) approach, a polymer is modeled 
as a sequence of N beads connected by elastic springs. The deformation of the bead-spring chain is then 
followed along the trajectory of its center of mass. In homogeneous and isotropic turbulence, Watanabe 
and Gotoh^ have shown that N = 2 beads are in fact sufficient to describe the stationary statistics of both 
the extension and the orientation of polymers, i.e., the deformation of a polymeric chain is dominated by 
its slowest oscillation mode. An analogous conclusion has been reached by Terrapon et alJ^ in a study of 
polymer dynamics in a turbulent channel flow. The model consisting of only N = 2 beads is known as 
the finitely extensible nonlinear elastic (FENE) dumbbell modeli The molecular approach is suitable for 
studying the deformation of passively transported polymers^r— . Two-way coupling molecular simulations 
are a recent achievement- 2 2 ~— and have not yet been employed in practical applications owing to their 
computational cost. In the molecular approach, indeed, the feedback of polymers on the velocity field is 
given by the forces exerted on the fluid by a huge number of polymeric chains 

For the above reason, in practical applications, numerical simulations of turbulent flows of polymer solu¬ 
tions generally use the continuum approach. The conformation of polymers is then described by means of 
a space- and time-dependent tensor field, which represents the average inertia tensor of polymers at a given 
time and position in the fluid. Such tensor is termed the polymer conformation tensor. An evolution equa¬ 
tion for the conformation tensor may in principle be derived from the FENE dumbbell model. Such equation, 
however, involves the average over thermal fluctuations of a nonlinear function of the polymer end-to-end 
vector; a closure approximation is therefore required. Peterlin^ proposed a mean-field closure according 
to which the average of the elastic force over thermal fluctuations is replaced by the value of the force at 
the mean-squared polymer extension. The resulting model was subsequently dubbed the FENE-P models 
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Within the FENE-P model, the back-reaction of polymers on the flow is described by a stress tensor field, 
which depends on the polymer conformation tensor. This continuum model is thus suitable for simulating 
turbulent flows of polymer solutions; it indeed amounts to simultaneously solving the evolution equation 
for the polymer conformation tensor and the Navier-Stokes equation with an additional elastic-stress term. 
The FENE-P model is widely employed in numerical simulations of turbulent drag reduction and has been 
successfully applied to channel flows,—shear flows , 31 i 32 and two- and three-dimensional homogeneous 
and isotropic turbulence— r— Nevertheless, although it qualitatively reproduces the main features of tur¬ 
bulent drag reduction, the FENE-P model generally does not yield results in quantitative agreement with 
experimental data (e.g. Ref. El). It is therefore essential to assess the validity of the assumptions on which 
the model is based. 

For laminar flows, the Peterlin approximation has been examined in detail (see Refs. El and El and 
references therein). In particular, the FENE-P model is a good approximation of the FENE model in steady 
flows, whereas appreciable differences appear in time-dependent flows. This observation suggests that in 
turbulent flows important differences between the two models should be expected. Several studies have 
subsequently investigated the validity of the Peterlin approximation in turbulent flows by comparing one-way 
coupling simulations of the FENE and FENE-P models— ~ 13 i 15 i 19 These studies have clearly shown potential 
differences among FENE and FENE-P models together with high sensitivity on the statistical ensemble and 
dependency on the degree of homogeneity of the underlying velocity fluctuations. We undertake a systematic 
analysis of the Peterlin approximation in three-dimensional homogeneous and isotropic turbulence by means 
of one-way coupling Lagrangian simulations of the FENE and FENE-P models. The size of our statistical 
ensemble (128 x 10 3 fluid trajectories and 2 x 10 3 realizations of thermal noise per trajectory) allows us 
to fully characterize the statistics of polymer extension and orientation. When the flow is turbulent two 
independent effects are at the origin of the discrepancies between the FENE and FENE-P models: one is 
directly related to the closure approximation for the elastic force, while the other is of a statistical nature 
and is a consequence of deriving the statistics of polymer deformation from that of the conformation tensor. 
By isolating these two effects, we compare the steady-state statistics and the temporal correlation of the 
extension and of the orientation of polymers in the FENE and in the FENE-P model. 

The rest of the paper is organized as follows. In Sect. [TT] we briefly review the FENE and FENE-P 
models. The Lagrangian simulations are described in Sect. IIIII The results of the simulations are presented 
in Sect. hyi Finally, conclusions are drawn in Sect. 


II. FENE AND FENE-P MODELS 


In the FENE model, a polymer is described as two beads connected by an elastic spring, i.e., as an elastic 
dumbbelli If the fluid is at rest, the polymer is in a coiled configuration because of entropic forces and its 
equilibrium extension is determined by the intensity of thermal fluctuations. If the polymer is introduced in 
a moving fluid and the velocity field changes over the size of the polymer, then the polymer can stretch and 
deform. The dynamics of the polymer thus results from the interplay between the stretching action of the 
velocity gradient and the elastic force, which tends to take the polymer back to its equilibrium configuration. 

The maximum extension of the dumbbell is assumed to be smaller than the Kolmogorov scale, so that 
the velocity held changes linearly in space at the scale of the dumbbell. The drag force on the beads is 
given by the Stokes law. Moreover, inertial effects and hydrodynamical interactions between the beads are 
disregarded. Polymer-polymer hydrodynamical interactions are also disregarded under the assumption that 
the polymer concentration is very low. Thus, the separation vector between the beads, R(t), satisfies the 
following stochastic ordinary differential equation (the FENE equation) 


dR 

dt 


a(t)R 



( 1 ) 


where R = |i?|, &ij{t) = djUiit) is the velocity gradient at the position of the center of mass of the dumbbell, 
£(t) is three-dimensional white noise, Rq is the polymer root-mean-square equilibrium extension, and t p is 
the polymer relaxation time ( t p is the time scale that describes the exponential relaxation of ( R 2 {t )) to 
its equilibrium value in the absence of flow). The three terms on the right-hand-side of Eq. ([1]) represent 
the stretching by the velocity gradient, the restoring elastic force, and thermal noise, respectively. The 
function / determines the elastic force and, in the FENE model, it has the following form: 


1 


/(C) 


1-C/T 2 ’ 


( 2 ) 
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where L is the maximum extension of the polymer. The elastic force diverges as R approaches L ; hence 
extensions greater than L are forbidden. Note that R(t) is a random vector and that, when u is turbulent, 
two independent sources of randomness influence its evolution: thermal noise and the velocity gradient itself. 

The polymer conformation tensor is defined as Cij = ( RiRj )j, where (-)j denotes an average over thermal 
fluctuations. To derive the evolution equation for C, we apply the Ito formula to RiRj and use Eq. ©: 


dR 


dRj i?.n 


Rl 


lift 


—(RiRj) — -jj-Rj+Ri — —f —Sij — aik(t)RkRj+ajk(t)RkRi-—f (R )RiRj+-^-Sij+^j [Rj£i(t) + Ri£j(t)]. 


dt 


dt 


(3) 

Whereas there is no Ito-Stratonovich ambiguity for Eq. ©, Eq. 0 should be understood in the Ito sense. 
We now average Eq. 0 with respect to the realizations of £(f) and make use of the following property of 
the Ito integral: (Ri^j(t))^ = 0. We thus obtain: 

— (RiRj ){ = aik(t)(RkRj)z + aj k (t)(R k Ri)t — — [( f(R )RiRj)( ~ Ro$ij\ ■ (4) 

Equation 0 is not closed with respect to C because of the term: 

Aj = (f(R 2 )RiRj)^. ( 5 ) 

To obtain a closed equation, Peterlin^ proposed the following approximation: 

Aj « f((R\)(RiRj\ = /(trC)Cy. (6) 

The resulting evolution equation for the polymer conformation tensor (the FENE-P equation) is: 

^ = a(t)C p +C p a\t) - 1 [/(tr C P )C P - R 2 0 l\ , (7) 


where I is the identity matrix and C p denotes the polymer conformation tensor calculated according to the 
Peterlin approximation. If the flow is turbulent, both C and C p have a random behavior. In the following, 
we shall denote: 


T=/(trC)C and T P = /(trC p )C P . (8) 

Equation 0 describes the evolution of the conformation tensor of a polymer along the Lagrangian trajectory 
of its center of mass. Numerical simulations of drag reduction use the Eulerian counterpart of Eq. 0, which 
is obtained by replacing dC p /dt with dtC p + u ■ VC P and a(t) with the Eulerian velocity gradient. In 
principle, the evolution equation for the conformation tensor should be coupled with the Navier-Stokes 
equations through an additional stress term proportional to V • T p Here, however, we focus on the impact 
of the Peterlin approximation upon the statistics of polymer deformation and consider passive polymers only 
(one-way coupling). In the rest of the paper, we thus study the relation between Eq. 0 and Eq. 0 when 
cr(t) is given by the incompressible Navier-Stokes equations in the turbulent regime. Some considerations 
are useful to guide our study: 

1. R(t) cannot be calculated from the solution of the FENE-P equation (Eq. 0 ). By contrast, C(t) can 
be constructed from the solution of the FENE equation (Eq. 0 ) by averaging the dyadic Ri(t)Rj(t) 
over the realizations of the noise £(i). The tensor C(t) can then be compared with C p (t ); 

2. because of the Peterlin approximation (Eq. 0), the FENE and FENE-P equations yield a different 
evolution for the conformation tensor. This holds for both laminar and turbulent flows; 

3. If the flow is turbulent, the statistics of R(t) differs from that of C(t), even if C(t ) is calculated from 
Eq. 0 (and hence no closure approximation is required). Consider for example the random variables 
R(t) and p(t) = yj ( R 2 (t))$ = y/trC(t). In general, the probability density function (PDF) of R(t) is 

different from that of p(t), as can be seen by noting that ((i?“)j) CT ^ ((R 2 )^)a- (a ^ 2), where (•)„• 
denotes an average over the statistics of the turbulent velocity gradient. 
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In conclusion, the FENE and FENE-P models differ for two reasons: the Peterlin approximation and the 
statistical effect due to the fact that the statistics of R(t) cannot be deduced from that of C(t). Hence, in 
the turbulent regime, the proper way to examine the Peterlin approximation is to first construct C (t) from 
the FENE equation and then compare its statistics with that of the solution of the FENE-P equation. If 
the statistics of R(t) is directly compared with that of C p (f), the error due to the Peterlin approximation is 
combined with the statistical effect discussed at point 3 above. This fact seems to have been overlooked in 
previous studies. 


III. LAGRANGIAN SIMULATIONS 

The dynamics of polymers is studied by using a database of Lagrangian trajectories that was previously 
generated to examine the dynamics of both tracer and inertial particles in turbulent flows * 40 : 41 The turbulent 
velocity field is obtained by direct numerical simulation of the three-dimensional incompressible Navier- 
Stokes equations: 

dtu + u • Vu = — Vp + z/V 2 w + /, V • u = 0, (9) 

where p is the pressure field and v = 2 x 10 -3 is the kinematic viscosity. The forcing f is such that the spectral 
content of the first low-wavenumber shells remains constant in time. The domain is a three-dimensional 
periodic box of linear size 2i r. Equations (0 are solved by means of a fully dealiased pseudospectral algorithm 
with second-order Adams-Bashforth time stepping. The number of grid points is 512 3 , while the integration 
time step is dt = 4 x 10~ 4 . In this simulation, the Kolmogorov time is t v = 4.7 x 10~ 2 and the Tavlor- 
microscale Reynolds number is R\ = 185 (for more details on the numerical simulation, see Refs. |40j and 
0). We expect our results not to depend significantly on the value of R\ except for some residual effects 
induced by intermittency in the statistics of the velocity gradients.— 

As mentioned in Sect. hd the inertia of polymers is negligible. Furthermore, their thermal diffusivity is 
very small compared to the turbulent diffusivity. Hence the center of mass of a polymer moves like a tracer 
and its position x c (t) satisfies the following equation: 

(j nr* 

~£ = u{x (10) 

Equation m is once again solved by using a second-order Adams-Bashforth scheme; a tri-linear interpola¬ 
tion algorithm is used to determine the value of the velocity field at the position of the polymer.— ^ After the 
statistically stationary state is reached for both the fluid motion and the translational dynamics of polymers, 
the positions of the center of mass of 128 x 10 3 polymers are recorded every At = 10 dt = 4x 10 -3 « t v / 10. 
The total integration time after steady-state is T = 13.2, which corresponds to 6 eddy turnover times 
approximately. 

The velocity gradient cr(i) is recorded every At along the trajectory of the center of mass of each polymer 
and is inserted into Eqs. 0 and 0 in order to determine the dynamics of the separation vector and 
of the conformation tensor. Equation © is solved by using the semi-implicit predictor-corrector method 
introduced by Ottinger^ the integration time step is equal to At for all values of the parameters. The 
initial condition for Eq. © is such that -R;(0) = Rq/VS, i = 1,2,3. Equation © for C p is integrated by 
means of the semi-implicit algorithm proposed in Ref. I13l which ensures that trC p ^ L 2 . 

The Weissenberg number is defined as Wi = t p /t. v and is the ratio of the time scales associated with 
the elastic force and with the velocity gradient. In our simulations, Wi varies between 10~ 2 and 10 2 . (An 
alternative definition of the Weissenberg number uses the maximum Lyapunov exponent of the flow, A, to 
estimate the reciprocal of the stretching time associated with the velocity gradient. In our simulations, 
A w 0.14Thus, the Weissenberg number based on the Lyapunov exponent is Wi* = Ar p « 0.14Wi). 
The squares of the equilibrium and maximum extensions of the polymer are Rq = 1 and L 2 = 3 x 10 3 , as in 
Refs. and[l9l. The number of realizations of thermal noise per Lagrangian trajectory is 2 x 10 3 . 

Finally, the statistics of polymer deformation is collected over the Lagrangian trajectories, over the real¬ 
izations of thermal noise, and over time (only for times greater than the time required for R(t) to reach the 
statistically steady state). 

We note that the statistics of the separation vector R in isotropic turbulence has been studied thoroughly 
by Watanabe and Gotoh.— The results on the statistics of R given below agree with those presented in 
Ref. 0 Here, we compare the statistics of C with that of C p , in order to determine the effect of the Peterlin 
approximation on the dynamics of polymers. 
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FIG. 1. Contour plot of P a (trA, trT) for Wi = 1 (top left panel), Wi = 7 (top right panel) and Wi = 20 (bottom 
left panel). The bottom right panel shows the averages of tr A (black squares) and tr T (blue circles) as a function 
of the Weissenberg number. 


IV. RESULTS 


In this Section, we examine the statistics of polymer extension and orientation in the FENE and FENE-P 
models. Before presenting the results, it is useful to define some notations. If the statistics of a random 
variable depends both on thermal noise and on the velocity gradient (as for instance in the case of R), 
its PDF is denoted as If the statistics of a random variable (e.g., trC) only depends on the velocity 
gradient, then its PDF is denoted as P a . The auto-correlation function of a scalar random variable X(t) is 
denoted as Fxit ) and the correlation time of X(t ) is: Tx = f 0 Fx{t)dt. The auto-correlation function of 
a statistically isotropic, unit random vector X(t ) is defined as: 


F x (t) = 2(\X(t + t 0 )-X(t 0 )\)-l (11) 

and the associated correlation time is defined as: Tx = / 0 °° Fx(t)dt. We denote by e, (* = 1,2,3) the 
unit eigenvectors of the rate-of-strain tensor S = (Vm +Vm t )/ 2; the unit eigenvectors e 7 ; are ordered by 
descending eigenvalue, i.e. e.\ is associated with the largest eigenvalue of S and e 3 with the smallest one. 
The direction of vorticity is u). Finally, we denote by Z\ and zf the first unit eigenvector of C and C p , 
respectively. 
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FIG. 2. Autocorrelation function of tr*4 (black squares) and trT (blue circles), for (a) Wi = 1, (b) Wi = 7, (c) 
Wi = 20. Panel (d) shows the correlation times T^a (black squares) and Thr (blue circles) rescaled by t v as a 
function of Wi. 


A. The Peterlin approximation 

As mentioned in Sect. ED the Peterlin approximation consists in replacing A with T in the evolution 
equation for the polymer conformation tensor (.4 and T have been defined in Eqs. (0 and (0 and, up to 
a concentration-dependent multiplicative constant, are the polymer stress tensors in the FENE and in the 
FENE-P models, respectively). The eigenvectors of A and T are the same, but their eigenvalues may differ. 
Indeed, Jensen’s inequality yields: tr„4 ^ trT. A first indication of the effect of the Peterlin approximation 
is thus given by the joint PDF P a (ti A, trT) (see Fig. [lj. For small polymer extensions or for small values 
of Wi, tr A and trT are approximately the same because 1 — R 2 /L 2 k, 1; hence the Peterlin approximation 
holds very well. By contrast, for large extensions or for large values of Wi, the deviations of trT from tr A 
are significant, i.e. the Peterlin approximation is poor. Notwithstanding, (tr T) a and (tr Tl)^- do not differ 
appreciably (see the bottom, right panel in Fig.QJ. This fact demonstrates that the study of average values 
may not suffice to investigate the validity of the Peterlin approximation. The full statistics of the separation 
vector and the conformation tensor should be investigated, which requires following the Lagrangian dynamics 
of a large number of polymers. 

In addition, the qualitative behavior of the temporal autocorrelations of tr A and tr T are similar (Figs. 0a) 
to0c)), but for large Wi the correlation time of tr A is shorter than that of trT (Fig.0d)). 


B. Statistics of polymer extension 

The statistical properties of the separation R are well understood. For small values of Wi, most of 
polymers are in the coiled state, i.e. their extension is close to the equilibrium one. Accordingly, the PDF of 
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FIG. 3. PDF of R/L (black squares), p/L (blue circles), p p /L (red triangles) for (a) Wi = 2, (b) Wi = 7, (c) 
Wi = 20. Panel (d) shows (R)^ rT /L (black squares), (p) a /L (blue circles), (p p ) CT /L (red triangles) as a function of 
the Weissenberg number. 




FIG. 4. PDF of R/L (black squares), p/L (blue circles), p p /L (red triangles) in semi-logarithnric scale for (a) Wi = 2, 
(b) Wi = 7, (c) Wi = 20. 


R has a pronounced peak at P 0 . As Wi increases, polymers unravel and become more and more extended. 
The transition from the coiled to the stretched state occurs when the Lyapunov exponent of the flow exceeds 
1/2 t p , i.e. at Wi> = l/24i (Note that some authors define the Weissenberg number in terms of the time 
scale associated with the exponential relaxation of sj ( R 2 (t )) instead of ( R 2 {t)) and hence obtain a critical 
Weissenberg number equal to unity.) 

At intermediate extensions, the PDF of R has a power-law behavior, i.e. P^ jCr (i?) ~ P _1_ “ for Rq -C 
R <C L 4^ This property of P^ i<T (i?) indicates that polymers with very different extensions coexist in the 
fluid; whether the coiled or the stretched state dominates depends on the value of Wi. The exponent a is 














































FIG. 5. Autocorrelation function of R (black squares), p (blue circles), p p (red triangles) for (a) Wi = 0.5, (b) 
Wi = 7, (c) Wi = 20. Panel (d) shows the correlation times Tr (black squares), T p (blue circles), T p (red triangles) 
rescaled by t v as a function of the Weissenberg number. 


positive in the coiled state and decreases as a function of Wi ^ 20 i 44 i 45 As long as a > 0, the FENE dumbbell 
model reaches a steady-state even as L tends to infinity (the L —> oo limit of the FENE model is known 
as the Hookean modcli). However, when a vanishes a steady-state PDF of the extension no longer exists 
if L —> oo. This behavior is interpreted as the coil-stretch transition^ Finally, if Wi increases beyond the 
value of the coil-stretch transition, a becomes negative and the maximum of Pj j(T (i?) moves from close to 
Rq to close to The statistics of R is shown in Figs. [3] and H] for different values of Wi. 

We noted in Sect. [TT] that the comparison between the FENE and the FENE-P model ought to be done in 
terms of the conformation tensors C and C p (rather than in terms of R and C p ). Let us denote p{t) = \JtvCit) 
and p p it) = y / trC p (t). To examine the influence of the Peterlin approximation on the statistics of polymer 
extension, we calculate p from the solution of Eq. © and p p from Eq. ©. We then compare P a (p/L) and 
P a (p p /L) in the steady state for different values of Wi. The plots shown in Figs. [3] and 0] correspond to the 
coiled state (Wi = 2), the coil-stretch transition (Wi = 7), and the stretched state (Wi = 20). 

For small values of Wi, P a (p/L) and P a (p p /L) do not differ significantly, because the extension of most 
of polymers is near to Rq and hence 1 — R 2 /L « 1 and C ~ C p (Fig. ©a)). For intermediate and large values 
of Wi, the statistics of R is characterized by a broad distribution of extensions around the mean value. 
The differences between P a (p/L) and P a (p p /L), therefore, are considerable. In particular, since (I? 2 ){ jCr is 
appreciably less than L 2 even for the largest Wi and /(C) rapidly diverges for C close to L [see Eq. ©], the 
restoring term in Eq. © is weaker than that in Eq. ©. Hence, in the FENE-P model large extensions are 
more probable than in the FENE model to the detriment of small and intermediate extensions (Fig. ©b) 
and (c)). As a consequence, the FENE-P model overestimates the average extension of polymers for large 
values of Wi (Fig. ©d)). An analogous behavior of the average extension has been observed in turbulent 
channel flowsi 13 ! 15 

We also note that whereas for large Wi the PDFs of R and p are approximately the same (Fig.©c)), they 
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FIG. 6 . PDF of alignment with (a) ei, (b) <= 2 , (c) e 3 and (d) oj for Wi = 20. 


are significantly different for intermediate or small Wi (Figs. [SJa) and (b)). Indeed, in the former case, the 
stretching action of the velocity gradient is very strong compared to the effect of thermal fluctuations, and 
in most realizations R ~ \JtrC. In the latter case, the effect of thermal fluctuations cannot be disregarded 
and the differences in the statistics of R and p (see Sect. d become evident. Furthermore, P^ jCr (I?/ L) 
and P a (p/L) mainly differ for small and intermediate extensions, because the large extensions are obtained 
in those realizations in which the velocity gradient is very intense and thermal noise can be neglected. 
The above results demonstrate that comparing the statistics of R directly with that of C p may lead to 
wrong conclusions; indeed, at small Wi, P^ jCT (P/i) and P a (p p /L) are clearly different, whereas P a (p/L) 
and P<j(p p /L) are close. 

The autocorrelation function of the extension is approximately exponential both in the FENE and in the 
FENE-P model. However, F p p ( t ) is a good approximation of F p (t) only for small Wi (Fig.jSja)). In addition, 
the FENE-P model captures the critical slowing down of polymers near the coil-stretch transition ) 17 ! 46 ! 47 
but for large Wi it underestimates the correlation time of the extension (Fig. [Sjd)). Once again, we note 
that, for small Wi, a direct comparison between P/j(f) and P p p(f) would lead to wrong conclusions about 
the effect of the Peterlin approximation on the temporal statistics of polymer extension. 


C. Statistics of polymer orientation 

In a coupled simulation of turbulent drag reduction, the feedback of polymers on the flow is of a tensorial 
nature; therefore, it depends not only on the extension of polymers but also on their orientation. 

For small values of Wi, the separation vector R exhibits a weak alignment with the eigendirections of the 
rate-of-strain tensor or with the vorticity^ because thermal fluctuations (which are isotropic) strongly affect 
the dynamics. In contrast, for intermediate or large values of Wi, thermal fluctuations have a negligible 
effect on the dynamics of polymers. Consequently, the evolution of the extension R decouples from that of 
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FIG. 7. Autocorrelation function of R (black squares), z\ (blue circles), if (red triangles) for (a) Wi = 0.5, (b) 
Wi = 7, (c) Wi = 20. Panel (d) shows the correlation times T^ (black squares), Tz 1 (blue circles), Ttp (red triangles) 
rescaled by t v as a function of the Weissenberg number. 


the orientation R = R/R 1 which behaves like the orientation vector of a rigid rod and is the solution of the 
Jeffery equation^ For Wi > 5, the statistics of R is therefore independent of Wi and coincides with that 
of the orientation of a rigid rod, i.e. R exhibits a moderate alignment with e 2 and a strong alignment with 
Cj (see Refs. [ft, Tf| and 143 and Fig. [6|). 

In the FENE-P model, the first eigenvector of C p , i.e. zf, can be interpreted as the orientation of the 
polymer, provided Wi is sufficiently large. The effect of the Peterlin approximation on polymer orientation 
can then be studied by comparing the statistics of z\ with that of the first eigenvector of C. Fig. [G] shows 
that the FENE-P model qualitatively reproduces the statistics of polymer orientation but underestimates 
the level of alignment. 

The autocorrelation function of R decays exponentially (Figs. EI a ) to[7[c)), which for large Wi is in 
agreement with the results for the autocorrelation of the orientation of a rod42 The correlation time of 
R increases for small values of Wi and saturates at large Wi (Fig. Eld)), because for large values of the 
Weissenberg number the dynamics of the orientation of R becomes independent of Wi. The FENE-P model 
agrees with the FENE model at small Wi, but quantitative discrepancies appear at large Wi (Figs. El a ) 
to Etc)). In particular, the FENE-P model underestimates the correlation time of polymer orientation 

(Figs. Efd). 


V. CONCLUSIONS 

Numerical simulations of turbulent flows of polymer solutions use the FENE-P model, which is based on 
the elastic dumbbell model but requires a closure approximation for the elastic term. We have examined 
the effect of the Peterlin approximation on the steady-state statistics of the extension and the orientation 
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of polymers. The FENE-P model captures the qualitative properties of the statistics, but for large Wi it 
underestimates the steady-state probability of small extensions and overestimates the probability of large 
extensions. As a consequence, the Peterlin approximation yields a greater average extension as well as a 
greater probability that polymers break under the action of a turbulent flow/ 0 To quantify this effect, one 
would need to couple the dynamics of polymers with a fragmentation model connected to the accumulated 
(or instantaneous) stress along each trajectory and to estimate the relative breaking rate.— Since large 
polymer extensions are more likely in the FENE-P model than in the FENE model, we also expect that the 
Peterlin approximation yields a stronger feedback of polymers on the flow in two-way coupling simulations 
of homogeneous isotropic turbulence with polymer additives. A similar argument, however, does not carry 
over to inhomogeneous flows like channel flows. In this case, indeed, drag reduction is caused by the strong 
stretching of polymers in the near wall region rather than by the dynamics in the bulk of the channel, where 
the flow is homogeneous and isotropic and a lower degree of stretching is observed ^di r 14 i 21 

As regards the temporal statistics of the end-to-end separation vector, both the correlation times of the 
extension and the orientation of polymers are underestimated by the FENE-P model. The FENE-P model 
also underestimates the level of alignment of polymers with the eigenvectors of the rate-of-strain tensor and 
with the direction of vorticity. 

It would be interesting to check to what extent these properties of the FENE-P model influence the 
dynamics of a polymer solution by comparing two-way coupling simulations of the FENE-P model in which 
the stress tensor is either calculated according to the Peterlin approximation or from molecular dynamics. 

In this paper, we examined the Peterlin approximation, as this is the main assumption in the construction 
of a continuum model of polymer solutions. However, it is worth recalling that the FENE-P model is based on 
the dumbbell model and hence on a very simplified coarse-grained description of a polymer macromolecule. 
Other approximations may thus impact the performance of the FENE-P model and its comparison with 
experiments. For instance, even for simple laminar flows, the dumbbell model reproduces the experimental 
observations only if the maximum extension and the effective bead radius are used as free parameters to fit 
the experimental data— ^ 
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